********************************************************************************
* Purpose: Longer run labor market outcomes where we have data until 2022. 
* Measured in quarterly outcomes. Regressions and coefficient plots
* For mothers and fathers (only mothers reported in paper)
********************************************************************************
do "D:\Data\workdata\708296\Project - Mother groups and mental health\Do\global.do" 

use "$work\final_data1", clear

* quarterly outcomes to simplify 
foreach v in barsel sick work barself sickf workf ui uif {
egen `v'_q1 = rowtotal(`v'1-`v'13)
egen `v'_q2 = rowtotal(`v'14-`v'26)
egen `v'_q3 = rowtotal(`v'27-`v'39)
egen `v'_q4 = rowtotal(`v'40-`v'52)
egen `v'_q5 = rowtotal(`v'53-`v'65)
egen `v'_q6 = rowtotal(`v'66-`v'78)
egen `v'_q7 = rowtotal(`v'79-`v'91)
egen `v'_q8 = rowtotal(`v'92-`v'104)
egen `v'_q9 = rowtotal(`v'105-`v'117)
egen `v'_q10 = rowtotal(`v'118-`v'130)
egen `v'_q11 = rowtotal(`v'131-`v'143)
egen `v'_q12 = rowtotal(`v'144-`v'156)
egen `v'_q13 = rowtotal(`v'157-`v'169)
egen `v'_q14 = rowtotal(`v'170-`v'182)
egen `v'_q15 = rowtotal(`v'183-`v'195)
egen `v'_q16 = rowtotal(`v'196-`v'208)
egen `v'_q17 = rowtotal(`v'209-`v'221)
egen `v'_q18 = rowtotal(`v'222-`v'234)
egen `v'_q19 = rowtotal(`v'235-`v'247)
egen `v'_q20 = rowtotal(`v'248-`v'259)
}
foreach v of varlist barsel_q13-barsel_q16 sick_q13-sick_q16 work_q13-work_q16 ui_q13-ui_q16 barself_q13-barself_q16 sickf_q13-sickf_q16 workf_q13-workf_q16 uif_q13-uif_q16 {
	replace `v' = . if mfr_yob > 2016
}	
foreach v of varlist barsel_q17-barsel_q20 sick_q17-sick_q20 work_q17-work_q20 ui_q17-ui_q20 barself_q17-barself_q20 sickf_q17-sickf_q20 workf_q17-workf_q20 uif_q17-uif_q20{
	replace `v' = . if mfr_yob > 2015
}
eststo clear
forvalue i = 1/20 {
eststo: reghdfe work_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\work_q`i'", replace)
}
eststo clear
forvalue i = 1/20 {
eststo: areg barsel_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\barsel_q`i'", replace)
}
eststo clear
forvalue i = 1/20 {
eststo: areg sick_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\sick_q`i'", replace)
}
eststo clear
forvalue i = 1/20 {
eststo: areg ui_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\ui_q`i'", replace)
}
eststo clear
forvalue i = 1/20 {
eststo: areg workf_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\workf_q`i'", replace)
}
eststo clear
forvalue i = 1/20 {
eststo: areg barself_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\barself_q`i'", replace)
}
eststo clear
forvalue i = 1/20 {
eststo: reghdfe sickf_q`i' peer_treat if tested_avg >= 1, absorb(time fe mfr_first) vce(cluster group) 
parmest, format(estimate min95 max95 %8.3f p %8.1e) saving("$work\sickf_q`i'", replace)
}

********************************************************************************
* Figure C8: Long run effects of exposure to a depressed peer on parental leave, 
* self-sufficiency, sick leave and unemployment benefits of mothers, quarterly status measure 
********************************************************************************
use "$work\ui_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\ui_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\ui1_20q.png", replace

use "$work\barsel_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\barsel_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

replace stderr = round(stderr,0.001)
twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\barsel1_20q.png", replace

use "$work\barself_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\barself_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\barself1_20q.png", replace

use "$work\work_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\work_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\work1_20q.png", replace

use "$work\workf_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\workf_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\workf_20q.png", replace

use "$work\sick_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\sick_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\sick1_20q.png", replace

use "$work\sickf_q1",clear
gen m = 1
forvalue v = 2/20 {
append using "$work\sickf_q`v'"
replace m = `v' if m == .
}
keep if parm == "peer_treat"

twoway (scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) /// 
(scatter estimate m, msize(small) msymbol(Oh) mcolor(red) alwidth(none) lcolor(black) lwidth(thick)) ///
(line estimate m, fcolor(black%10) alwidth(none) lcolor(red) lwidth(thick)) ///
(line max95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)) ///
(line min95 m, fcolor(black%30) alwidth(none) lcolor(gs10) lwidth(medium) lpattern(dash)), $gsq20 legend(off) xtitle("Quarters after birth") ytitle("Difference") ylab(, tlcolor(gs10) labgap(6-pt) format(%4.2f)) yline(0, lcolor(gs10) lwidth(medthick))
graph export "$out\sickf_20q.png", replace

